The window at the edge of chaos in a simple model of gene interaction networks 
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As a model for gene and protein interactions we study a set for molecular catalytic reactions. 
The model is based on experimentally motivated interaction network topologies, and is designed to 
capture some key statistics of gene expression statistics. We impose a non-linearity to the system 
by a boundary condition which guarantees non-negative concentrations of chemical concentrations 
and study the system stability quantified by maximum Lyapunov exponents. We find that the 
non-negativity constraint leads to a drastic infiation of those regions in parameter space where the 
Lyapunov exponent exactly vanishes. We explain the finding as a self-organized critical phenomenon. 
The robustness of this finding with respect to different network topologies and the role of intrinsic 
molecular- and external noise is discussed. We argue that systems with inflated 'edges of chaos' 
could be much more easily favored by natural selection than systems where the Lyapunov exponent 
vanishes only on a parameter set of measure zero. 

PACS numbers: 87.16.Yc, 64.60.Ht, 82.39.Rt, 89.75.Da, 



I. INTRODUCTION 

Most complex systems, living systems in particular, 
are characterized by remarkable degrees of stability and 
at the same time by a tremendous potential of flexibil- 
ity and adaptability. This has led some authors to de- 
fine cornplex and living systems as living at the 'edge of 
chaos', [l|, [2, y, 3 meaning -in a somewhat picturesque 
way- that it takes only tiny changes in the system to 
move it from a stable and regular mode into a chaotic 
phase where large portions of phase space can get sam- 
pled. The concept is that systems at the edge of chaos 
are especially well suited for adaptation and information 
processing in the sense that adaptability is associated by 
the possibility of finding adequate new states in possibly 
changed environments at very fast rates. It has been ar- 
gued that living systems at the edge of chaos would get 
favored by natural selection, and that life has evolved to- 
wards such a special region in parameter space [3|. In 
many dynamical systems the edge of chaos is a very spe- 
cial set of points in parameter space, often of measure 
zero, characterized by the system's maximal Lyapunov 
exponent A passing through zero. It is not clear how 
systems can get regulated towards (or have evolved to- 
wards) such a limited set of critical points, even though 
some interesting ideas have been proposed in this direc- 
tion 0. Even in the simplest maps like the logistic map, 
the dynamics exactly at these special points can become 
highly non-trivial [6|. 

It is evident that living systems have evolved towards 
stable systems in stationary disequilibrium. Various au- 
thors argue that a key principle of living systems is their 
ability to replicate [7]; corresponding rate equations for 
molecular replicators have been proposed for a long time, 
beginning with [g] . As such, basic molecular reactions in 
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living systems (e.g. protein production or degradation) 
have to be autocatalytic. If autocatalytic reactions are 
not balanced by degradation and/or thermostatic net- 
flow of substance to and fro the system (like in a flow 
reactor), concentrations of molecular products will di- 
verge in the replicator. A stationary state can be estab- 
lished when production and decay (flow) rates of inter 
cellular molecules effectively balance each other, p. [lO|. 
In this sense stability (stationarity) provides a natural 
selection criterion. Catalytic reactions are simply de- 
scribed by reaction networks, which contain production 
and degradation rates. Given current developments in 
genomics- and proteomics technology some facts about 
these networks become known. By now there is some 
evidence that these (directed) networks show scale-free 
(SF) topological organization Tl|, Il3l. On the basis of a 
given molecular reaction topolog y 11311 several gene net- 
work models have been proposed [ij, llSl . Il6l | . In principle 
two different approaches have been pursued: discrete ap- 
proaches, using Boolean networks [l7| and continuous ap- 
proaches, using ordinary or stochastic differential equa- 
tions [l^, [l^, I20I. l2li . Combinations of both have also 
been reported [22, [23| . 

Recently the importance of noise in molecular reaction 
networks has been stressed and its relevance has been ex- 
perimentally demonstrated [2J, [23, [2a| • For example the 
level of noise can determine whether cells in Drosophila 
become epidermal or neural cells [23] • Further it was 
shown that low reproduction rates of DNA and impor- 
tant regulatory molecules forbid to neglect stochastic ef- 
fects |i28i] . Intrinsic noise, microscopic events within the 
cell, and extrinsic noise, such as cell to cell variations, 
are experimentally distinguishable [29| . In this context a 
stochastic differential equation model has been proposed 
for regulatory transcription networks [30| ■ 

In this work we study a simple linear, noise driven 
dissipativc model for catalytic molecular reactions. We 
impose a non-linearity to this system by assuming that 
concentrations can not become negative. We demon- 



strate that this non-hnearity changes the 'the edge of 
chaos' from a point where A = 0, to extended regions of 
vanishing Lyapunov exponents. The model offers a way 
to understand how systems naturally evolve and adapt 
towards a widened 'edge of chaos'. 



II. THE MODEL 

We assume that gene to gene interactions can be mod- 
eled as chemical reactions between proteins, mRNA and 
other nucleic material. Let us denote the concentrations 
of proteins i at time t by the vector Pi{t) and of RNA 
molecules j by rj (t) . For convenience let us combine all 
types of concentrations into a single vector x = (j>,r). 
The size of the vector (number of products) we denote 
by N. The simplest linear model to capture all possible 
interactions is given by 
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where A^ is the matrix of reaction rates. Even though 
this model is clearly an over-simplification of reality it has 
been frequently used recently [15, jil, 32]. Let us assume 
that these rates are not perfect constants but fluctuate 
according to A^^ (t) — Aij + ^ij (t) , for example through 
thermal noise. For simplicity let ^y- be an iid process 
with zero mean. Replacing A'^ by A we get 
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Regardless of the distribution of t,ij , and assuming that 
X will converge to a reasonably stationary distribution, 
according to the central limit theorem, the sum of the 
right hand side will yield a random number from a Gaus- 
sian distribution, which we shall denote by r\i G iV(0, a). 
For simplicity we shall assume S^a also Gaussian, i.e. 
^ii = £,i € N{0, (T)with the same variance a Vz. We need a 
final addition of our model Eq. ^ , to incorporate a fur- 
ther experimental observation. Many gene-products (e.g. 
mRNA levels) fluctuate around some characteristic value, 
x^, across the cell cycle. In Fig. [T^ we show the expres- 
sion levels for mRNA levels of 10 randomly picked genes 
from the yeast genome (S.cerevisae) over 2 cell cycles at 
17 time points taken at 10 minute intervals |33l. |34| . The 
number of measured genes in the budding yeast genome 
was N = 6220. To incorporate these characteristic levels 
we chose values a;^ from some distribution. In Fig[T|3 we 
show the experimental distribution of mRNA expression 
levels of the yeast genome, defined as the time-average 
over cell cycles, xf = {xi{t))t- In the following x" is 
taken from uniform, Gaussian or the experimental distri- 
bution. This fixes our model to be 



dt 
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A%{xj-x^j)+i,{t)x, + T], 



(3) 



with ^i E N{0, a) and rji G iV(0, a) the multiplicative and 
additive noise components respectively. Multiplicative 
and additive noise can be interpreted as intrinsic and 
extrinsic noise as used e.g. in [29|. Both intrinsic and 
extrinsic noise are present in gene networks. 

To be able to interpret x as concentrations we have to 
introduce the constraint 



x,{t) >0 \/i 



(4) 



which means that regardless of Xi of ([3]), Xi(t) can never 
be below zero. This imposes a non-linearity onto the 
system and makes it non-trivial. 



A. The interaction matrix 

Before solving the system we have to specify the inter- 
action network, i.e. the matrix elements (chemical rates) 
of A. It is obvious that the network will be directed and 
weighted. Diagonal elements An < are decay rates 
of the products, off-diagonal rates Aij can be both pos- 
itive or negative, corresponding to activation or inhibi- 
tion. Further, it is clear that first not all products can 
interact with each other, i.e. a large number of matrix 
elements will be zero, and second most rates are not avail- 
able from experiments. We are thus led to model A as a 
random matrix in the following way. Using terminology 
from network theory the 'degree' ki of product i is de- 
fined as the number of products that can be regulated by 
product i. The class of interaction networks can now be 
specified by the 'degree distribution'. There is evidence 
that protein networks [i2| and metabolic networks [1J| 
are scale-free networks p{k) ^ k~'' , characterized by a 
degree distribution with average degree (k) > 4, and an 
exponent 7 ~ 2.2. In the following we generate such net- 
works and contrast them to classical random networks, 
i.e. Erdos-Renyi graphs (ER) [35l ? ] with the same av- 
erage degree. If the number of non-zero rates in A is 
denoted by L, the average connectivity is (fc) = L/N. 

Once it is decided which products interact with each 
other, i.e. Aij ^ 0, the actual rates have to be fixed. 
We assume these being Gaussian, Aij G iV(0, cr^). This 
is supported experimentally e.g. by [32], where a least 
squares fit of synthetic gene network models to real 
data indicates that the normal distribution of interac- 
tion weights gives the best results. The (negative) decay 
rates An we take constant and equal Vi in the following. 



B. A note on multiplicative noise 

Note that the diagonal component of Eq. Q - ignoring 
the positivity constraint ~ immediately reminds at the 
stochastic differential equation. 



dt 



x = f{x)+g{x)at)+v{t) 



(5) 



with / = —AiiX and g{x) = x. This Langevin process 
has been exactly solved [3a|, the solution (probability 
distribution of x) being a g-exponential, 



p{x) ^\l + {q~ l)(ix^ 



li/(i-«) 



(6) 



with /3 — {—An + (T/2)/cr, where (1 — q)^^ is the asymp- 
totic power exponent. In Fig. Ic we show experimen- 
tal data confirming the power law aspect of niRNA con- 
centrations in yeast in the same data [33| . l^Xi{t) = 
Xi{t) — a;i(t — 1) is the difference in gene expression lev- 
els between two consecutive measurements; P{> Ax) is 
the cumulative distribution, for all i and t. In the same 
plot we show results of a numerical simulation of the 
model Eq. Q with the N = 6000 ER topology for A 
at {k) — 20, for the two cases: first, a — Q and ct > 
(Gaussian noise model), and second a > and a — Q 
(multiplicative noise model). Further a g-exponential fit 
[42] to the data is shown (broken line), with an effective 

Qc.uin -^ 1.00 . 



C. stability 

Biological systems are sufficiently stable, i.e. products 
are not produced ad infinitum, and at the same time are 
sufficiently dynamical. If our model is reasonable we have 
to find the regions in parameter-space where such non- 
trivial stability is ensured. If we ignore for a moment the 
positivity condition, Eq. ([4]), and the stochastic terms 
in Eq. ([3]), the stability of the system is dominated by 
the largest real part of the eigenvalues of the interaction 
matrix A. If there are no non-negative real parts of the 
eigenvalues, the system will be asymptotically stable. If 
the distribution of off-diagonal elements in A is normal 
[33 | with variance a\, the eigenvalue spectrum is - ac- 
cording to a powerful result from random matrix theory 
- a circle in the complex plane (Girko's circular law), 



For a fully connected matrix, with L 



N' 



zero entries in A the radius of this circle p is equal to 
the product of the standard deviation and square root 
of the system size N . For non-fuUy connected networks, 
L < iV^, the radius is given by (see e.g. ^881] ) 



CTA 



VL/N = aAV{k) 



(7) 



If the diagonal elements of the random matrix are from 
a zero-mean distribution, Girko's circle is centered at the 
origin of the complex plane. In our case we have An < 
and the center of the circle will be shifted to the position 
(-A,„0), seee.g. [H. 

Now, including the positivity constraint and the 
stochastic dynamics, it is obvious that the eigenvalue 
spectrum of A will only be a part of the story. To define a 
measure for system stability that captures these aspects, 
maybe the simplest choice is the maximal Lyapunov ex- 
ponent 
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FIG. 1: (a) A set of 10 randomly picked gene-expression 
trajectories of yeast {S. cerevisiae) over two cell cycles [331 ]. 
Stationary state values x" are defined as the time-averages of 
gene expression levels, xq levels corresponding to the shown 
genes are indicated by horizontal lines, (b) Distribution of 
stationary state values a;". Inset shows the cumulative distri- 
bution, (c) Cumulative distribution of gene expressions incre- 
ments for the same data (circles) for the numerical simulation 
of the evolution of gene expression data with two evolution 
models with additive Gaussian (boxes) and multiplicative (di- 
amonds) external noise. 



where Sx{t) = x{t) — x'{t), is the difference between two 
trajectories, where x'{t) results from a small perturbation 
in the initial condition, ||(Sx(0)|| ^ 1. 

For the system without the positivity condition the A 
can now be related to p in the following way 



(8) 



Xm)-Pm)~Au = aAV{k)-A, 



(9) 



where the An is the eigenvalue spectrum shift discussed 
above. 

For the case of the full model i.e. with the positivity 
condition, we hypothesize the following scenario: With 
strong noise levels after some time several trajectories 
diffuse to hit zero. For long enough times the expected 
number of these trajectories will be half of the total num- 
ber, N/2. This amounts to a reduction of system size by 
one half (with fixed connectedness), i.e. N -^ Ncs = N/2 
and L -^ L^s = L/A. For connectivity this means, 
(fc) -^ (fccff) — {k)/2. We would thus expect the asymp- 
totic (large (fc)) behavior of A of the full model as a func- 
tion of connectivity 
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For small (fc) where no (or few) trajectories hit zero, of 
course, we expect Eq. (j9|) to hold. More generally, it 
is reasonable to assume that for given connectivity and 
noise levels, there will be an effective connectivity. 



(fceff) = {k)Non/N 



(11) 



where Non is the number of trajectories not at zero. 

Let us finish with commenting on a potential stabi- 
lizing role of multiplicative noise [40|. Consider the one 
dimensional case of our model 



—-x — a(x ~ Xo) + ^x + ri 
at 



(12) 



with ^ e N{0,a) and rj G N{0,a). The evolution of a 
perturbation Sx thus follows 



-—Sx 
dt 



aSx + S.Sx 



with the solution 



(13) 



(14) 



The Lyapunov exponent is proportional to a—a^/2 show- 
ing that the system can be stable even for positive a. 



III. RESULTS 

We numerically solve the model Eq. ([3]), now with 
positivity condition and compare with the above predic- 
tions. For numerical simulations we generated scale- free 
(SF) and ER networks of sizes of iV = 200, 500, 1000. To 
vary (k) we adjusted the number of non-zero rates L in 
the matrix A. For scale-free networks we fixed the scal- 
ing exponent 7 — 2.2. All the following results are aver- 
ages over 20 random realizations of networks for a given 
parameter set. The Lyapunov exponents were computed 
from datasets of 1000 timepoints, after discarding the ini- 
tial 200 timesteps. xq was chosen from the experimental 
distribution of Fig. [TJd. We do not observe noteworthy 



FIG. 2: Maximum Lyapunov exponents A as a function of 
average degree (k), averaged over 20 realizations, for ER net- 
works. The influence of the positivity condition on forming 
a plateau is immediately seen. Simulations are shown with 
(a = a — 0.1) and without {a — a = 0) noise. A'' = 500, 
An = -0.4. Lines are Eqs. JS} and (fT0|) . 



changes of results when using unifornr or Gaussian dis- 
tributions. In Fig. [21 we show the solution for A for the 
ER network, as a function of (k), with (triangles) and 
without (circles) positivity condition. The correspond- 
ing theoretical predictions Eqs. ^ and (fTO|) . are drawn 
as broken and solid lines, respectively. The case with- 
out positivity condition is completely explained by the- 
ory over the entire range of {k), Eq. ([9]). In the case 
with the constraint the asymptote follows Eq. (|10p . as 
expected. It is also seen that for small (fc), Eq. ([9]) is 
valid. The main finding of this paper is that within a 
(fc)-window between about 10 and 30 a plateau forms 
where A practically vanishes. The constraint dynamics 
allows a full range of connectivities to support a 'life at 
the edge of chaos'. 

The stability of the system for different network 
topologies, sizes and various noise components is shown 
in Fig. [3l Figure [3K indicates that both, network size 
and degree distribution are slightly influencing the width 
of the plateau. While in the (fc) -^ N region there is no 
signiflcant difference in system stability, the low connec- 
tivity region shows a size effect on the A = plateau. 
The effect of network topology is relatively small, the 
curve pertaining to SF always being slightly below the 
ER networks, see Fig. [3^. While the width of the 
plateau is always wider for the random distribution of 
links, in the (k) — > region, the system is more sta- 
ble (smaller A) for SF networks. For higher connectivity 
regions ((fc) > 30) the difference between random and 
scale-free networks becomes indistinguishable due to nu- 
merical inaccuracy. Figure [31d shows the influence of pure 
multiplicative (a > 0, cr — 0) and pure additive noise 
(ct = 0, cr > 0) on the plateau width, compared to the 
deterministic process, a = cr = 0. With multiplicative 
noise the plateau becomes significantly wider, while ad- 
ditive noise hardly shows any effect when compared to 
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FIG. 3: Lyapunov exponents for same parameters as in 
previous figure, for (a) different network types and sizes 
A'' = 200, 500, 1000, (b) noise effects for multiplicative noise 
{a = O.I and ct = 0), additive noise {a — and a = O.I), 
compared to the deterministic process (ct = a = 0). (c) A 
compared to the number of inactive nodes as a function of con- 
nectivity. ER, iV = 1000, An = -0.4, and cr = CT = O.I. (d) 
Clearly <ja is not a constant, and declines with inactive nodes, 
as expected. ER, iV = 1000, An = -0.4, and a = ct = 0. 



the deterministic process. Plateau widths are collected 
in Tab. 1. 

To understand the formation of the plateau (A = 0), 
one could naively expect Eq. © to hold at the plateau 
with the modification that (k) is replaced by (/ccff) from 
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TABLE I: Zero-A plateau widths A for an ER network, with 
A'' = 500. The width A is defined as the region of connectivity 
where |A| < 0.005. For the situation of variable multiplicative 
noise the additive noise was fixed to o = O.I, for the variable 
additive noise, a = 0.1. Cases for different An are shown. 



Eq. pTjl . This then would imply 



iVon = 



NAl 



\{k) 



(15) 



In Fig. [5}; we show that the tail of A^on ~ — (fc) and does 
not follow the naive expectation from Eq. (fT5|) . ^ {k)~^ ■ 
The key to understanding the formation of the plateau 
lies in the fact that aA is not a constant. When A ^ 0, 
and nodes start becoming inactive, this amounts to effec- 
tively reducing the interaction matrix. For each product i 
hitting zero, this means that row and column of matrix A 
can be dropped. (In biological terms, once the concentra- 
tion of some product type reaches zero, these molecules 
stop playing a role in the regulation of other products, 
the effective regulation network consisting just of active 
nodes {xi > 0), gets smaller.) The key observation is that 
this does not leave the variance of the matrix elements 
Aij unchanged, but systematically reduces the variance 
the more nodes become inactive. This can be seen as a 
selection mechanism in which the most active reactions 
(largest reaction rates in A) will hit the boundary first, 
and will be the first ones to be removed on average. This 
mechanism drives the system to a self-organized critical 
state at A = 0. We show the situation in Fig. [3Jl; (ta 
declines with the number of inactive nodes N — Nq^, as 
expected. With further increasing (fc) the number of in- 
active products saturates at half the network size and Eq. 
(Uni) holds. 



IV. DISCUSSION 

We have studied the stability of a simple stochastic 
model of catalytic reaction equations for cellular prod- 
ucts such as mRNA molecules or proteins. The system is 
driven by intrinsic molecular noise (multiplicative) and 
external (additive) noise. We show that the model cap- 
tures some essential experimental features, such as the 
fat tail distribution of concentration changes. Imposing 
an intuitively natural constraint on the system, (non- 
negativity of concentrations) we observe the forming of a 



plateau of vanishing Lyapunov exponents. The dynami- 
cal stability of concentrations in catalytic regulatory net- 
works, defined with Eq. ([3]) has three extended phases in 
parameter space (here connectivity). In the first phase 
the system is asymptotically stable, A is negative. After 
being exposed to a random perturbation in this phase 
the system always relaxes to its steady state x^. The 
main finding of this work is the existence of a second 
phase, where A ~ extends over extended regions. This 
is in marked contrast to the dynamics of many other 
non-linear systems, which have A = only at a set of 
points. The emergence of this phase can be fully ex- 
plained within the model. At higher connectivities some 
products will start to reach the boundary. Those prod- 
ucts with the largest reaction rates will - on average - 
hit the boundary first. These are then removed from the 
system. This means that the variance of the effective 
reaction rates will get driven downward as a function of 
connectivity. The variance of rates and the number of 
active nodes balance each other to exactly arrive at the 
critical point A = 0. This is a self-organized critical ef- 
fect. A physical analogy is the evaporation at boiling 



temperature, where molecules of higher-than-average en- 
ergy leave the liquid first, keeping the (critical) temper- 
ature at the boiling point. The third phase is defined by 
A > where system is dynamically unstable and concen- 
tration levels are diverging. We studied the dependence 
of the A = plateau on two network topologies, ER and 
SF, where a remarkably small effect was found. We found 
that with multiplicative noise the size of the plateau can 
be varied while additive noise showed relatively little ef- 
fect. For strong enough noise levels of either type the 
plateau breaks down. 

In [4l| it was noted that neural networks can perform 
most complex computations if the dynamics of random 
threshold gate networks is at the critical boundary be- 
tween the ordered and chaotic regime. If we interpret 
gene-regulatory networks as computing devices perform- 
ing hundreds of optimization problems simultaneously, it 
is plausible that evolution would have selected among the 
most efficient variations - working at the edge of chaos. 

This project was supported by WWTF Life Science 
Grant LS 139, and by the Austrian Science Fund FWF, 
project P19132. 
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